Draft version September 18, 2012 

Preprint typeset using I^T^X style cmulatcapj v. 12/16/11 



RHAPSODY: I. STRUCTURAL PROPERTIES AND FORMATION HISTORY FROM A STATISTICAL SAMPLE 

OF RE-SIMULATED CLUSTER-SIZE HALOS 

Hao-Yi Wu, 1,2 Oliver Hahn, 1 Risa H. Wechsler, 1 Yao-Yuan Mao, 1 Peter S. Behroozi 1 

iRavli Institute for Particle Astrophysics and Cosmology; Physics Department, Stanford University, Stanford, CA 94305 
SLAC National Accelerator Laboratory, Menlo Park, CA 94025 
2 Physics Department, The University of Michigan, Ann Arbor, MI 48109; hywu9umich.edu 
Draft version September 18, 2012 

ABSTRACT 

We present the first results from the Rhapsody cluster re-simulation project: a sample of 96 "zoom-in" 
simulations of dark matter halos of lO 14 8±O O5 ft. _1 M0, selected from a 1 ft~ 3 Gpc 3 volume. This 
simulation suite is the first to resolve this many halos with ~ 5 x 10 6 particles per halo in the 
cluster-mass regime, allowing us to statistically characterize the distribution of and correlation between 
halo properties at fixed mass. We focus on the properties of the main halos and how they are affected 
by formation history, which we track back to z = 12, over five decades in mass. We give particular 
attention to the impact of the formation history on the density profiles of the halos. We find that the 
deviations from the Navarro-Frenk- White (NFW) model and the Einasto model depend on formation 
time. Late-forming halos tend to have considerable deviations from both models, partly due to the 
presence of massive subhalos, while early-forming halos deviate less but still significantly from the 
NFW model and are better described by the Einasto model. We find that the halo shapes depend 
only moderately on formation time. Departure from spherical symmetry impacts the density profiles 
through the anisotropic distribution of massive subhalos. Further evidence of the impact of subhalos is 
provided by analyzing the phase-space structure. A detailed analysis of the properties of the subhalo 
population in Rhapsody is presented in a companion paper. 

Keywords: cosmology: theory — dark matter — galaxies: clusters: general — galaxies: halos - 
methods: N-body simulations 



1. INTRODUCTION 



Galaxy clusters are powerful probes of cosmological pa- 
rameters and have played a key role in the development 



of th e current ACDM paradigm (see, e.g., Allen et al. 



2011 for a review). For example, the spatial distribution 
and abundance of galaxy clusters reflect the growth rate 
of large-scale structure and the expansion rate of the 



energy (e.; 



Vikhlinin et al. 2009 



Rozo et al.|2010afc, neutrino mass 



Mantz et al.| 


2010b! 


., Mantz et al 


|2010a| 



and the validity ot g eneral relativity on 
cosmic scales ]&!g., |Rapetti et al.|2010 2012). In the near 
future, the massive influx ot mufti-wavelength data (e.g., 
SPTj^ACTpJ PlanclQ eRositfFJ PanSTARRSpJ DES*| 
Euchcj^J LSST|^]) will greatly enhance the sample size of 
galaxy clusters and reduce the statistical uncertainties 
in cluster cosmology. However, the constraining power 
of galaxy clusters will depend on how well various sys- 
tematic uncertainties can be controlled, including the 
relations betwee n observable properties and mass (e.g., 
Rozo et aL||2012 ); the robustness of cluster identification 
and centering (e.g., Rykoff et al.||2012 ); and the effect of 



1 The South Pole Telescope; http://pole.uchicago.edu/ 

2 Atacama Cosmology Telescope; http:/ /www. princeton.edu/act/ 

3 http://www.esa.int/planck 

4 Extended ROentgen Survey with an Imaging Telescope Array; 
http:/ /www. mpe.mpg.de/eROSITA 

5 The Panoramic Survey Telescope & Rapid Response System; 
http:/ /pan-starrs. ifa.hawaii.edu/ 

6 The Dark Energy Survey; http://www.darkenergysurvey.org/ 

7 http://sci.esa.int/euclid/ 

8 The Large Synoptic Survey Telescope; http://www.lsst.org/ 



White et al.|2010| ) 



viewing angle and projection (e.g., 

One essential way to understand these systematic uncer- 
tainties is through N-body simulations, which have been 
appli ed to study galaxy clusters for more than a decade 
(e.g. |Tormen et al.|1997| jMoore et al ||1998l [Gh igna et al. 
1998[ also see Kravtsov fc B organi 2 012| tor a more gen- 



erai review). In the era of large-sky survey and precision 
cosmology, it is desirable to have controlled simulation 
samples that can help us understand the statistical dis- 
tribution of the properties of galaxy clusters and the 
correlation between observables, as well as their detailed 
structures and evolution. Since massive galaxy clusters 
are rare, cosmological simulations need to cover a large 
volume to include a fair number of these sy stems (e.g., the 
MultiDark simulation |Prada et a l.||2011| a nd the recent 

3c ■ 



Millennium XXL simulation Angulo et al.||2012| ) . How- 



ever, given limited computational resources, the detailed 
substructures of halos are not well-resolved in these sim- 
ulations. Instead of using a cosmological volume with a 
single resolution, one can focus on particular systems and 
re-simulate them with higher resolution. This so-called 
"zoom-in" technique provides a powerful way to study indi- 
vidual chistersysternsin detail in a cosmological context 



(e.g.,|Tormen et al.|1997 


Moore et al.|1999||Navarro et al. 


2004fCao et al. 2U05 


bleed et al. 2005). Nevertheless, 



small number of cluster-size syste ms (e.g., the cu rrent 
high-resolution Phoenix simulation Gao et al.|2012 ) and 
galactic halos (e .g., the Via Lactea 11 simulat ion Die 
Imand et al.||2008 and the Aquarius simulations |Springei 



|et al.||2008| ). Therefore, few statements have been made 
about the statistical properties of well-resolved subhalos 
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in the mass regime of galaxy clusters. 

In this work, we perform re-simulations of a large num- 
ber of cluster-forming regions in a cosmological volume 
(side length 1 h~ 1 Gpc) to create a high-resolution statis- 
tical cluster sample, Rhapsody, which stands for "Re- 
simulated HAlo Population for Statistical Observable- 
mass Distribution studY" . The current sample includes 
96 halos of mass 10 14 ' 8 /i _1 M0 with mass resolution 
1.3 x 10 8 /i -1 M Q . One of the main goals of Rhapsody 
is to create a sample of cluster-size halos at fixed mass 
that enables us to make statistical statements about the 
halo population that is relevant for current and imminent 
cluster surveys. 

In Figure 111 we compare Rhapsody 8K (main sample) 
and Rhapsody 4K (a factor of 8 lower in mass reso- 
lution) to sever al N-body simulations in th e literature 
(Mille nnium II |Boylan-Kolchin et al .||2009|; Millenium 
XXL [Angulo et al.||2012 



MultiDark IPrada et al 



Bolshoi Klypin et al. 2011 



2011 ; Consuelo and Carmen 



from LasDamas; McBride et al. in preparation]; P hoenix 



Gao et al.||2012| ; Aquarius |Springel et al.||2008| ). Ha- 
os from zoom-in simulations are presented by symbols, 
while halo populations inside cosmological volumes are 
presented by curves with the shape of the halo mass func- 
tion. Rhapsody is in a unique regime in terms of the 
number of halos in a narrow mass bin simulated with high 
particle number. It is also worth noting that Rhapsody 
is currently the largest sample of halos with more than 
a few times 10 6 particles per halo at any given mass in 
the literature. Our repeated implementation of the re- 
simulation method makes the simulation suite statistically 
interesting and computationally feasible. 

The Rhapsody sample is highly relevant to several 
current observational programs. Fo r example, the galaxy 



cluster catalogs from the SDSS (Koester et al. 2007 
Hao et al. 2010) include many tens ot thousands ot 



photometrically-selected clusters and have provi ded a 



rich sample for multi- wavelengt h mass calibration (Rozo 
et al.|2 010b Rykoff ct al. 20121, cosmological constraints 



^ozo et al. 2007, 2010a), and studies of the cluster galaxy 
populations ( |Hansen et al.||2009[ ). The Cluster Lensing 
And Supernova su rvey with Hubble Mu lti-Cycle Treasury 
Program (CLASH; Postman et al. 2011 1 focuses on 25 mas- 
sive clusters and aims to establish unbiased measurements 
of cluster mass-concentration relation of these clusters. 



Recently, von der Linden et al. (20121 published accurate 
weak lensing mass calibrations of 51 massive clusters, fo- 
cusing on understanding various systematics for cluster 
count experiments. In addition, various X-ray programs 
have been efficiently identifying massive cluster s; for ex- 



ample, the ROSAT Brightest Cluster Sample (Ebeling 
let al.|2000|), the ROSA T-ESO Flux-Limited X-ray sample 
jBohringer et al.||2004| ), and the MAssive Cluster Survey 
( jEbeling et al.||2010[ ). These samples have achieved high 
completeness and provided a relatively unbiased selection. 
Relatively recently, massive galaxy clusters have also been 
detec ted through the Suny aev-Z erdovich (SZ) effect by 
ACT jMarriage et aL|2Q]T|, SPT (|Williamson et al.|201l| ), 
and Planck ( |Planck Collaboration et al.||2011[ ). These 
detections have ushered in an era ot high-purity detection 
of high-redshift galaxy clusters. The Rhapsody sample 
is in a mass regime similar to these observational pro- 
grams and can provide a statistical description of the 
dark matter halos associated with these clusters. 
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Figure 1. Comparison of the halo samples in various N-body 
simulations; Rhapsody is in a unique statistical regime of well- 
rcsolvcd massive halos. The number of halos (per 0.1 dex in mass) 
is shown as a function of number of particles inside the viriai 
radius of the halo. Symbols correspond to halos in re-simulation 
projects; the Rhapsody 4K and 8K samples are shown as two 
colored stars (M v i r = 10 14 ' 8±0 f)5 /i — ^^Mq). Curves correspond to 
halos in different cosmological volumes, and black stars on these 
curves correspond to the number of halos of the same mass as 
Rhapsody. We note that Consuelo and Carmen both include 50 
volumes, and only one volume is presented here. 

This paper presents the first results from the Rhap- 
sody simulations. We first characterize the formation 
history and the density profiles of the 96 main halos. We 
then explore how formation history impacts halo concen- 
tration and the deviation from the Navarro-Frenk- White 
(NFW) profile. We find that the deviations from the 
NFW model systematically depend on formation time 
and are impacted by the presence of massive subhalos. 

We connect the density profile to the phase-space struc- 
tures of halos and find that late-forming halos tend to 
have outflows within i? v ir, which can be also be attributed 
to massive subhalos. We have also investigated the shape 
parameters for the spatial distributions and velocities of 
dark matter particles. We find that the shape parameters, 
after removing massive subhalos, have no strong correla- 
tion with formation time, indicating that the deviation 
from spherical symmetry alone cannot account for the 
trend of profile and formation time. 

This paper is organized as follows. In <j2j we detail 
the simulations . We present the formati on h istor y of the 
main halos in |3.1| and merger rate in |3.2| In §4.1[ we 
present the density profiles and compare various fitting 
functions; in |4.2| we demonstrate the effect of formation 
history on the density profile. In f|5j we analyze the 
impact of formation time on the phase-space structure. 
The shapes of spatial distributions and velocities of dark 
matter particles and their alignments are discussed in fj6j 
We conclude in |7j 



In a second paper in this series ( Wu et al. 2012 hereafter 
Paper II), we will present the properties of subhalos in 
our sample and the impact of formation time on them, 
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which is more complex in the sense that the subhalo 
properties depend on the selection method of subhalos, 
the stripping experienced by a subhalo, and the resolution 
of the simulation. 

2. THE SIMULATIONS 
The Rhapsody sample includes 96 cluster-size halos of 



mass M,, 



10 



14.8±0.05 



h Mq, re-simulated from a cos- 



mological volume of 1 /i _1 Gpc. Each halo was simulated 
at two resolutions: 1.3 x 1O 8 /i _1 M (equivalent to 8192 3 
particles in this volume), which we refer to as "Rhap- 
sody 8K" or simply "Rhapsody"; and 1.0 x 1O 9 /i _1 M 
(equivalent to 4096 3 particles in this volume), which we 
refer to as "Rhapsody 4K." These two sets allow de- 
tailed studies of the impact of resolution. The simulation 
parameters are summarized in Table [l] 

The initial conditions were generated with the multi- 
scale initial condition generator Music (Hahn & Abel 



2011 ). The particles w ere then evolve d using the public 
version of Gadget-2 ( |Springel|[2005l ). The halo finding 



was p erformed with the pha se-space halo finder Rock 
STAR (jBehroozi et al.||2011a). Finally, merger trees were 



constructed with the 



jravitationally-consistent code of 



|Behroozi et al. (2011b). We provide more details on our 
methods below. 

All simulations in this work are based on a ACDM 
cosmology with density parameters ft m = 0.25, Q\ — 
0.75, fif, = 0.04, spectral index n s = 1, normalization 
er 8 = 0.8, and Hubble parameter h — 0.7. 

Figure [2] shows the images of 90 halos at z — in 
the Rhapsody 8K sample. Halos are sorted by their 
concentration and subhalo number, as described in the 
following sections. Figure [3] shows the evolution of 
4 individual halos, selected as extremes in the distri- 
bution of concentration and subhalo number. Movies 
and images for each individual halo are available at 
http : //risa. Stanford . edu/rhapsody/. 

2.1. The cosmological volume 

Our re-simulations are based on one of the Carmen 
simulations from the LArge Suite of DArk MAtter Sim- 
ulations (LasDamas; McBride et al.). The simulation 
represents a cosmological volume of 1 h~ 1 Gpc with 1120 3 
particle s. Its initial conditi ons are generated with the 
code of Crocce et al. ( 2006 1 based on the second-order 
Lagrangian perturbation theory (2LPT), and the N-body 
simulation was run with the Gadget-2 code. Rhapsody 
uses the same cosmological parameters as Carmen. 

When selecting targets for re-simulation from the mas- 
sive end of the halo mass function, we choose a mass bin 
that is narrow enough so that mass trends of halo prop- 
erties are negligible but at the same time wide enough 
to include a sufficient number of halos for statistical 
analyses. Here we focus on a 0.1 dex bin surrounding 



logioM v 



14.8. This mass range allows us to select 



100 halos in a narrow mass range, and is well-matched 
to the masses of the massive clusters studied in X-ray, 
SZ, and optical cluster surveys. 

2.2. Initial conditions 

For each of the halos in our sample, we generate multi- 
resolution ini tial conditions using the Music code (Hahn 
fc Abel|20TH ). We use the same white noise field of UAR- 
MEN (1024' 5 of its 1120 3 modes) to generate the large-scale 



perturbations consistent with Carmen. The equivalent 
resolution ranges from 256 3 in the lowest resolution region 
to 8192 3 (4096 3 for the 4K sample) in the highest resolu- 
tion region. In between, the mass resolution changes by 
factors of 8 every 8 times the mean inter-particle distance. 

For each of our re-simulation targets, we choose a zoom- 
in volume that is 40% larger than the Lagrangian volume 
of the friends-of-friends halo at z — 0. This choice has 
been tested to provide a well-converged dark matter den- 
sity profile in our convergence tests. With this setting, 
no low-resolution particle was found within the virial 
radius of any targeted halo. The typical number of high 
resolution particles per simulation is thus about 42/5.4 
million for 8K/4K with a standard deviation of 18%. 

In Music, particle displacements and velocities have 
been computed from the multi-scale density field using 
2LPT at a starting redshift of 49, in accordance with 
Carmen. The use of 2LPT is important for statistical 
studies of such massive systems since their m asses depend 
on the accuracy of the in i tial conditions (e. g., Crocce et al. 



2006 



2012 



Tinker et all[2008l [Reed et al.||2012| jBehroozi et"aT 



and McBride et al., in preparation) . 



2.3. Gravitational evolution 

After generating the initial conditions, we evolve each 
cluster-forming region using t he public version of the 
Gadget-2 code ( Springel|2005 1 . Gravitational forces are 
computed using two levels of particle-mesh together with 
the force tree to achieve a force resolution of comoving 
3.3/6.7 h~ 1 kpc in the Rhapsody 8K/4K for particles in 
the high resolution region. For each simulation, we save 
200/100 snapshots logarithmically spaced in scale factor 
a between a = 0.075 and a = 1 for the 8K/4K sample. 

We note that the virial masses of the re-simulated halos 
change somewhat with the improved resolution. As a 
result, a fraction of the halos fall outside the narrow tar- 
geted mass range logioM v ; r = 14.8 ± 0.05. In most cases, 
the masses scatter slightly upwards. We discard those 
halos falling outside the targeted mass bin of Rhapsody 
to keep the mass selection clean. In principle, to obtain 
all halos in the 14.8 ± 0.05 mass bin in the re-simulated 
sample, one needs to re-simulate a wider range of masses 
around 14.8 to include all halos that end up in the tar- 
geted bin. However, the large suite of re-simulations thus 
required is beyond the scope of this work. Thus, we note 
that Rhapsody does not strictly include the complete 
sample of halos within logioM vir = 14.8 ± 0.05 in either 
the original volume or the re-simulations. However, we 
do not expect this fact to affect the results presented in 
this paper, because the main approach in this paper is 
stacking all halos in Rhapsody for sufficient statistics 
and our sample should be unbiased. Global statistics for 
halos in this bin in the entire cosmological volume (for 
example, the two-point correlation function) are not used 
in the present work. 

2.4. Halo and subhalo identification 

Our simulations are post-processed with the adap- 
tive ph ase-space halo finder Rockstar ( jBehroozi et al 
2011a), which can achieve high completeness in finding 
subhalos e ven in dense environments ( see also Knebe 
et al.|[20lT ). Based on the phase-space information, small 
subhalos passing through the dense central region of the 
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Type 


Name 


Mass Resolution 

[/l^Mg] 


Force Resolution 
[fc-^kpc] 


Number of Particles 
in Simulation 


Number of Particles 
in Each Targeted Halo 


Full Volume 


Carmen 


4.94xlO IU 


25 


1120 a 


12K 


Zoom- in 


Rhapsody 4K 
Rhapsody 8K 


l.OxlO 9 
1.3xl0 8 


6.7 
3.3 


5.4M a /4096 3 (equiv.) 
42M a / 8192 3 1 (equiv.) 


0.63M b 
4.9M 6 



"The mean number of high-resolution particles in each zoom-in region. 

'The mean number of high-resolution particles within the R v j r of each targeted halo. 



Table 1 

Simulation parameters. 



main halo can be robustly identified. This feature is es- 
pecially important for studying the subhalo populations, 
which we focus on in Paper II. We note that the algo- 
rithm is only applied to high-resolution particles in the 
simulations. 

Rockstar pays special attention to major merger 
events (two halos of similar mass merge with each other), 
which arise frequently in the formation history of Rhap- 
sody halos (because of their high masses) and often cause 
difficulties in the construction of merger trees. During 
a major merger between two halos, a large fraction of 
dark matter particles appear as unbound to either of 
the merging halos, even though they are bound to the 
entire merging system. Therefore, regular unbinding pro- 
cedures tend to result in ambiguities or inconsistencies 
in halo mass assignment. Rockstar addresses this issue 
by computing the gravitational potential of the entire 
merging system, thus making the mass evolution of halos 
self-consistent across time steps. 

2.5. Merger trees 
We apply the gravita tionally-consisten t merger tree al 



gorithm developed by Behroozi et al. (2011b) to the 
200/100 output snapshots which were saved between 
z = 12.3 and z = for Rhapsody 8K/4K. The idea 
behind this new merger tree implementation is that the 
stochasticity in N-body simulations often leads to failures 
in halo finding. For example, the halo finder might find 
a spurious halo that is in reality a random density fluc- 
tuation at a certain time step, or the halo finder might 
miss a halo because it falls below the detection threshold 
at that particular time step. Given these limitations in 
halo finders, previous implementations of merger trees 
often encounter problems in linking halos across different 
time steps. The gravitationally-consistent merger tree 
algorithm resolves this issue by comparing adjacent time 
steps to recover missing subhalos and remove spurious 
halos, thereby improving the completeness and purity of 
the halo catalogs and ensuring correct linking of halos 
across time steps. This algorithm compares two adjacent 
time steps and can be summarized as follows: (1) It takes 
the halos at the later time step and evolves their positions 
and velocities backward in time, deciding whether the pro- 
genitors are missing or incorrectly linked. (2) It takes the 
halos at the earlier time step and looks for its descendant 
in the later time step. If the descendant is missing, the 
algorithm decides whether a merger occurs or the current 
halo is spurious. For details of the imple mentation, we 
refer the reader to Behroozi et al. (2011b). 



2.6. Summary of halo properties 

In Table [2] we summarize the key halo properties dis- 
cussed throughout this paper. In Figure [4j we present the 
distributions of and correlations between several of these 



properties. We use the rank correlations throughout this 
work to avoid the impact of outliers. 

In this work, the halo mass definition is based on the 
spherical overdensity of virialization, A vir , with respect 
to the critical density, /0 C rit- We use the center of the 
phase-space density peak calculated by Rockstar as 
the center of a halo. Based on this center, we draw a 
sphere with radius i? v ir so that the mean overdensity 
enclosed is equal to A v i r p cr i t . With the cosmological 
parameters used herein, A v , r = 94 with respect t o the 
critical density at z = ( Bryan fc Norman||1998[ ); i.e., 
A V i r = Ag4 C =A37@ m . The subscripts c and m indicate 
the overdensities with respect to the critical density and 
mean matter density. For reference, in Table [2] we list 
halo masses and radii based on several commonly-used 
overdensity values: A 20 0m = A 50 c, A 20 o c , A 500m = A 125c , 
and A 500c . 

In Table [2] we list two properties that are closely re- 
lated to halo mass: the maximum circular velocity and 
the velocity dispersion of dark matter particles. The max- 

that 



imum circular velocities is defined at a radius r max 
maximizes y/GM(< r)/r: 



<GM{< r n 



(1) 



The velocity dispersion is calculated based on dark matter 
particles: 



1 



y l — l 



(2) 



We note that the correlation between M v ; r and a v is 0.32, 
indicating that there is a non-negligible residual mass- 
velocity dispersion scaling despite the narrow mass range 
of our sample. For this reason, in the remainder of the 
paper, we have always verified carefully that the quoted 
correlations between various properties are not simply 
driven by mass. 

3. THE BUILDUP OF CLUSTER-SIZE HALOS 

In this section, we present the mass accretion history 
and merger rate of the main halos in Rhapsody, paving 
the way for further discussions of the impact of formation 
time on halo properties. The formation history of dark 
matter halos is known to correlate with various halo 



and subhalos (e.g., 


-n r>: — 

Wechsler et al.||2002| 


Zhao et al.|2003| 


Harker et al.| 


2006 


Hahn et al. 2007; Maulbctsch et al. 


2007 


Li et al. 


2008 


I. On the other hand, the merger rate 
os serves as a baseline for modeling 


ot dark matter ha 



several processes in galaxy formation, including the build- 
up of stellar mass and supermassive black holes, the star 
formation rate, the color and morphology evolution (e.g., 
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Figure 2. Images of 90 RHAPSODY halos at 2 = 0. The halos are first sorted by concentration (high concentration on the upper rows, low 
on the bottom). In each row, the halos are then sorted by the number of subhalos (selected with v max > 100 km/s, high number of subhalos 
on the left columns, low on the right). Each image has a physical extent of 4 h~ x Mpc on a side, which is slightly larger than the average 
virial radius of 1.8 7i -1 Mpc. 
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Figure 3. Evolution of four Rhapsody halos. From top to bottom, the images show the progenitors of four halos at z = 3 and 2 = 1, and 
the halo at z = 0. The four halos chosen are the corners of Fig 2. From left to right, they have high concentration, high subhalo number 
[337]; low concentration, high subhalo number [377]; high concentration, low subhalo number [572]; low concentration, low subhalo number 
[653]. Halo 572 has the highest concentration, the least late-time accretion, and the most dominant central halo of our full sample. It is also 
the halo with the most massive progenitor at z = 3. Each panel has a comoving extent of 4 /i -1 Mpc on a side, and is centered on the most 
massive progenitor in each case. 
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Figure 4. Distributions of and correlations between main halo properties and formation history parameters. The number in red in each 
panel shows the rank correlation coefficient, and its font size reflects the magnitude of correlation. 
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Property 


Median Frac. Scatter 


Mean Frac. SD 


Def. 


Af vir [ft- 1 *!©] 

fivir [/l _1 Mpc] 

(j v [km / s] 

Vmax [km/s] 


6.4 x 10 14 0.072 
1.8 0.024 
1,400 0.041 
1,300 0.036 


6.4 x 10 14 0.061 
1.8 0.028 
1,400 0.039 
1,300 0.042 




^200m 
^500m 

M 500c 


7.3 x 10 14 0.083 
5.9 x 10 14 0.065 
5.0 x 10 14 0.063 

3.4 x 10 14 0.13 


7.3 x 10 14 0.077 
5.9 x 10 14 0.056 
5.0 x 10 14 0.062 

3.4 x 10 14 0.12 




2.6 




^200m 
-R200c 

i?500m 
-R5OOC 


2.3 0.028 
1.6 0.022 
1.3 0.021 
0.84 0.043 


2.3 0.026 
1.6 0.019 
1.3 0.021 
0.83 0.04 


h-iMpc; \ 


2.6 




2l/2 

%>OL 

7-/3 


1.2 1.1 
0.58 0.44 
0.67 0.058 

2.3 0.22 


1.6 1.1 
0.56 0.43 
0.67 0.059 

2.2 0.23 






r s [/i _1 Mpc] 
C NFW 
CNFW-likc 

CEinasto 
7NFW-likc 
7Einasto 
^Einasto 


0.34 0.31 

5.3 0.26 
5 0.27 

4.9 0.23 

3.4 0.43 
3 0.13 

0.24 0.42 


0.36 0.27 
5.3 0.23 
5.1 0.32 
4.9 0.28 
3.9 0.34 
3 0.14 
0.27 0.41 






b/a 
c/a 
T 


0.75 0.12 
0.63 0.12 
0.71 0.27 


0.76 0.12 
0.63 0.12 
0.69 0.27 








0.82 0.096 
0.72 0.098 
0.17 0.32 


0.82 0.093 
0.72 0.094 
0.18 0.3 





Table 2 

Properties of Rhapsody halos at z = 0. The second column shows the median, the third column corresponds to the ratio of the 68% scatter 
to the median. The fourth column shows the sample mean and the fifth column corresponds to the ratio of the standard deviation to the 

mean. 
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Figure 5. Mass accretion history of the main halos in Rhapsody. 
The gray curves indicate the trajectories of individual halos; the 
black curve indicates the average over all halos. The cyan curve 
indicates the pseudo-evolution of a static halo. We show the virial 
mass of the most massive progenitor at each output redshift. The 
bottom panel shows the residual of the fit for three analytical 
models. 



Kauffmann & Hachnelt 2000; Hopkins et aLl|2006l [De 



Lucia fc Biaizot||2007[ |Behroozi et al.||2012p . In addition 



there has been an ongoing effort to meas ure merger rates 
in observat ions (e.g., Bell et al.|2006 Lotz et al.|2008 Xu 
et al. ||2012" ). While a study ot the implications tor galaxy 



formation is beyond the scope of this paper, we note that 
the merger rate provided here can be applied to modeling 
galaxy formation in massive clusters. 

3.1. Mass accretion history 

Figure [5] presents the mass evolution of the main halos 
in Rhapsody. For each main halo identified at z = 0, we 
search through its merger tree to find the most massive 
progenitor at each redshift. The gray curves show the 
evolution of M v ; r for individual halos, and the black curve 
shows the average of all halos. We note that the fractional 
dispersion is roughly constant for all redshifts. In the 
upper panel, we add a cyan line showing the pseudo 
evolution of a static halo (with a non-evolving density 
profile and M vir = 1O 14 - 8 /i _1 M at z = 0) expected 
simpl y from the evolution of A v j r and the critical density 
(e.g. jDiemer et al.|2012 |. 

Our halo formation histories span almost 5 orders of 
magnitude in mass and cover the 13.8 Gyr between z = 12 
and z — 0, allowing us to test the validit y of parameter- 
izations proposed in the literature (e.g., Wcchslc r et al. 
20021 jvan den Bos ch||2002l |Tasitsiomi et al P004[ |7haol 
et al.||2009| |McBride et al. |2009p over much larger ranges 
in both time and mass. We fit trie following three models: 



1. An exponential model ( Wechsler et al.|[2002 ) 

M(z) = M e- az . (3) 
We note that this model assumes a constant mass ac- 



cretion rate represented by the exponential growth 
index 

dlnM , , 

-^- =a - (4) 

The curves in Figure [5] do not follow straight lines, 
indicating a systematic deviation from the pure 
exponential model. For this model, a formation 
time proxy can be defined as 



ln2/a 



(5) 



2. An expon ential-plus-power law model with two pa- 
rameters ( |McBride et~aT|2009| ) 

(6) 



M(z) 

We note that 
dlnM 



dz 



M (l + zf 



7 — j3 when z « 1 



(7) 



Thus, 7-/3 can be used as a measure for the late- 
time accretion rate. Analogous to Equation [5j a 
formation time proxy can be defined as 



z 7 _ = In 2/(7-/3) 



(8) 



One can alternatively solve M(z^ 7 ) = Mo/2 numer- 
ically. However, we note that z^ 7 , z~—p, and 7 — /? 
are completely correlated with each other. 

3. An exponential-plus-power law model with one pa- 
rameter, motivated by Equation [6} When fitting for 
Equation [61 we observe that the two parameters /3 
and 7 are highly correlated (see Appendix [A} . We 
thus adopt a 1-parameter model that incorporates 
this correlation 



M(z)=M (l + z) 



-4.61+5. 277 -fz 



(9) 



The two numerical coefficients are obtained by an 
optimization scheme described in Appendix |A"] 

Our fitting procedure and the goodness-of-fit are also 
discussed in Appendix |X] The bottom panel of Figure [5] 
shows the average of the residual for the individual fit, 
(Mfit/-Mt ruc ). The pure exponential model does not pro- 
vide the curvature needed to fit the data. The two 
exponential-plus-power law models work almost equally 
well, but significant residual remains for z < 2. 

3.2. Merger rate 

The merger rate experienced by a dark matter halo as a 
function of time and merg er ratio is a specific prediction 
of the ACDM model (e.g., |Press fc Schechter|1974|[Lacey 
fc Cole|1993 ) and has been calibrated with ever- improving 
precision using the exten ded Press-Schechter model (e.g., 
Neistein fc Dekel 2008[) and N-body simulations (e.g., 
Fakhouri fc Ma||2(Mijl |Stewart et al.||2009| |Genel et al.| 
2009 ). In this work, we use the new m erger tree im- 
plementation of Behroozi et al.| 
the merger rate of cluster-size h, 



(2011b|) to re-examine 
alos. 



We consider the 
epoch of a merger event to be the earliest redshift where 
a smaller halo enters the virial radius of a larger halo, 
and the merger ratio, fi = M mel ^/M main , is defined at 
this epoch. 

Figure M presents the merger rate for Rhapsody ha- 
los. The left panel shows the cumulative merger rate: 
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2006 



Figure 6. The merger rate of RHAPSODY halos. Left: Cumulative merger rate as a function of rcdshift (per halo). The average number of 
merger events each main halo has experienced since a given z is shown for three different merger mass ratios fi = M mergms /M mam . Right: 
Differential merger rate as a function of merger mass ratio. The number of merger events per halo per dfi per dz is shown as a function of fi. 
The three different curves correspond to different redshifts, and the trend is independent of redshift. For both panels, each colored band 
corresponds to the standard deviation of the curve of the same color. 

The x-axis specifics a look-back redshift, and the y-axis 
shows the average number of merger events each main 
halo has experienced since this redshift. Different curves 
correspond to different merger mass ratio thresholds, y, > 
0.01, 0.03, 0.1, and 0.3, and the shaded regions indicate 
the standard deviation of the sample. We identify the 
redshift of the last major merger, Zi mm , using \i > 0.3. 
Our merger rate is con sistent with the t rend o f the main 
halo mass presented in 



an ongoing effort (e.g. , [Bullock et al.|2 001 Wc chsler et a l. 
INeto et al.|2007[|lVlacci6 et al.|2008| |Gao et al.|200 



Fakhouri et al. (20101, who have 
ations 1 and 11 and do not have 



used Millennium Simu 
sufficient statistics in our mass regime. 

The right panel shows the merger rate as a function of 
the merger mass ratio fj,. We plot the differential number 
of merger events each main halo has experienced, per d[i 
per dz, for a given merger ratio. The different curves 
represent different redshifts. We find that the merger 
rate trends are almost invariant with redshift. The blue 
shaded region corresponds to the scatter of the blue curve 
(z = 0.05) and indicates the large variation of merger 
rate from halo to halo. The logarithmic slope, -1.9 2, is 
very close to the valu e -2 in Fakhouri & Ma (|2008|) and 



Prada et al. 2011 Bhattacharya et al. 2011) and is of 
increasing importance tor interpreting observations. As 
mentioned in the introduction, the CLASH project is a 
major effort of the Hubble Space Telescope and aims for 
detailed and unbiased measurements of the density profile 
of galaxy clusters, which are tests of both the ACDM 
paradigm and our understanding of the assembly of clus- 
ters. In addition, the modeling of the concentration-mass 
relation imp acts the interpretat ion of the weak lensing 
results (e.g., King fc Mead||2011[) and X-ray results (e.g., 
Ettoriet al.||201o| . 



Fakhouri et al. 



( 2010 1. We also plot the fitting formula in 
Eaknouri et ai. (20101 as a black dashed curve (z — 0.05). 
While the overall slope agrees well, their normalization is 
slightly lower but still agrees within our error bar. 

Figure HI shows that Zi mm correlates with Z1/2, as well as 
various halo structural properties, which will be detailed 
in the following sections. 

4. THE IMPACT OF FORMATION TIME ON THE DENSITY 

PROFILE 

The density profiles of dark matter halos have been 
shown to foll ow the universal N avarro-Frenk-White 
(NFW) form ( Navarro et al.p"997 ) and can be well char- 
acterized by the concentration parameter c. Calibrating 
the concentration-mass relation and its scatter has been 



In this section, we first provide fits to the density profiles 
of the halos in the Rhapsody sample. The halos of galaxy 
clusters, being the most massive objects in the universe, 
assemble very late in cosmic history, so that recent violent 
merger and accretion events are expected to impact their 
density profiles. For this reason, we study in detail the 
relation between their formation history and profiles in 
the remainder of the section. 



4.1. Fitting the density profile 

The left panel of Figure [7] shows the density profiles 
of the main halos in Rhapsody. The black curve corre- 
sponds to the mean density profile, while the gray curves 
correspond to individual halos. We use 32 bins between 
10~ 2 - 5 i?vir and i?vir equally spaced in logr for plotting 
these curves. We note that the difference in the binned 
density profile between each re-simulated halo and its 
low-resolution version is between 5% and 10%. 

For each main halo, we use all dark matter particles 
between 13 /i~ 1 kpc and i? v ir to fit the density profile, 
adopting the maximum-likelihood estimation without bin- 
ning in radius (explained in detail in the Appendix |B[) . 
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Figure 7. Left: Density profiles of the main halos in RHAPSODY (gray: individual halos; black: average). The bottom panel shows the 
residuals with respect to three models (fits to the mean profile). Right: The impact of formation time on density profile. Halos are binned 
into quartiles based on their Z\ and the mean of halos in each quartile is shown in the upper panel. The bottom panel shows the residuals 
with respect to the NFW profile. Deviations from the NFW profile are systematic, with the largest deviation seen for late-forming halos and 
monotonous decrease when considering halos in the earlier forming bins. 



There has been an ongoing effort to quantify the deviation 
from NFW and to search for alternative fitting functions 
e.g. |Navarro et al.||2004l IMerritt et al.|2006 



we define 



Gao et al. 



2008. JNavarro et ai. 2010). In this work, we fit three 
parameterizations: 



R v 



^Einasto " 

7Eina S to = 24^0 (slope at R vir ) , 



r-2 

OBinasto 



(15) 
(16) 



The NFW profile 
p{r) 



{r/r s )(l + r/r s ) 2 
l + 3(r/r.) 



Pcrit 

dlxip 
dlnr 1 + (r/r s ) 



(10) 



(11) 



which is characterized by one parameter cnfw = 

Rvir I '8 • 

An "NFW-like" profile with a free outer slope 7 
p(r) $c 



where 7Einasto is the slope of the logarithmic density 
profile at i? v ; r and can be compared with 7NFW-iikc- 

The best fit values of these parameters are summarized 
in Table [2j and the comparison of the goodness-of-fit 
is detailedm Ap pendix [B} Fo r the N FW fit, our mean 
value agrees with jPrada et al" (2011) and Bhattacharya 
Our standard deviations are a(c)/c 



et al. 



Pcrit (r/r s )(l + r/r s )i- 
dlri p 1 + j(r/r s ) 
dhir 1 + (r/r s ) 



(12) 



(13) 



which reduces to the NFW profile when 7 = 3. 
This profile can be characterized by two parameters, 

7NFW-likc and CNFW-likc = (-RvirAs)(7NFW-like ~ 

2). The latter is defined so that i? v ir/cNFW-iiko 
equals the radius at which the density slope is —2. 
We impose r s < i? v ; r in the fitting procedure to 
avoid possible divergence of 7. 

3. The Einasto profile ( |Einasto|[l965l ) 

d\np = _ 2 /_r 

dlnr \ r -2, 

This model is characterized by two parameters, r_2 
and QfEinasto- To compare with the other two models, 



(14) 



|2011|. Our standard deviations are c(c)/c = 0.26 
and a (log j o c) = 0.11, which are slightly small er than the 
values quoted in Bhattacharya et al. ( 2011[ ) (0.33 and 
0.16 based on A 2 oo c ). Tins difference is presumably due 
the decreasing scatter in the concentr ation-mass relation 
with i ncreasing mass; as mentioned in Bhattacharya et al. 
(2011), their scatter is slightly smaller for massive halos 
(for~M 20 o > 8 x 10 14 r^o, the scatter is o(c)jc = 0.28, 
which is very close to our value). 

In the left panel of Figure [7J we add a bottom panel to 
compare the residual of these three models, (p) /pat- Here 
Pfa(r) is obtained by fitting the stacked binned density 
profile of all halos, and the legend shows the best-fit 
parameters for each model. We note that these values are 
slightly different from the average values of halos (shown 
in Table [2]); i.e., the fit of the average results in slightly 
higher concentrations than the average of the fit to each 
halo. Among these three profiles, the Einasto profile fits 
to the stacked density profile best, deviating by up to 
5%, whereas the NFW profile deviates by up to 10%. A 
similar trend of devi ations has also b een shown in the 
Phoenix simulations ( Gao et al.||2012[ ) and the Aquarius 
simulations ( Navarro et al.|[20TU[ r 
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4.2. Density profile and formation history 

We explore several aspects of the impact of the for- 
mation history on halo density profiles. In §4.2. 1| we 
present the correlation between concentration and forma- 
tion time, as well as the i mpact of formation time on the 
deviation from NFW. In [4.2.2 we show that the slope of 



profile also depends on formation time and that subhalos 
alter the slope profile of late-f orming halos. We show the 
evolution of concentration in j ]4.2.3| 

4.2.1. Formation time, concentration, and the deviation from 

NFW 

Figure [4] includes cnfw and the goodness-of-fit of NFW, 
defined as 



A N fw = max|M(< r) - M NFW (< r)\/M v 



(17) 



We find that both quantities are strongly correlated with 
z i/2, Zlmmi and 7-/3, and the correlation is strongest 
zi/2- The correlation between concentration and 



with 

formation time is we ll known (e.g., Navarro et al 



1997 



Wechsler et al.|[2002 ), and the canonical explanation is 



that the concentration of halos reflects the physical density 
of the universe at their formation epoch. The correlation 
between the deviation from NFW and formation time 
can be understood through the relaxedness of halos: late- 
forming halos tend to show larger deviation from NFW 
because of recent mergers and late-time mass accretion. 
Although highly concentrated clusters tend to be closer 
to NFW, we note that cnfw and Anfw only have a weak 
correlation. 

To further explore the correlation between formation 
time and concentration, as well as the deviation from 
NFW, we split our halos into quartiles based on Zi/%. 
In the right panel of Figure [7J we present the average 
density profile of halos of each quartile. The upper panel 
shows that the formation time can lead to systematic 
differences in the density profile. For each average density 
profile, we fit the NFW profile and present the residual, 
(p( r ))//°NFW- in the bottom panel. Deviations from 
NFW are again systematic; around 0.1 i? v j r NFW fits 
tend to overestimate. The latest-forming quartile tends 
to have the largest deviations from NFW, while the two 
early-forming quartiles show less deviation. We note that 
the deviation from NFW cannot be entirely attributed to 
the departure from spherical symmetry. As we will show 
in SjHJ the shape and triaxiality of halos do not strongly 
correlate with formation time and cannot account for the 
observed deviation from NFW. 

4.2.2. The slope of density profile: Impact of subhalos 

To further understand the impact of formation time 
on the deviation from NFW, we compare the logarithmic 
slope of the density profile in the quartiles of highest 
and lowest Zi/ 2 ■ The left panel of Figure JHl shows the 
absolute value of the logarithmic slope ofthe density 
profile, r = — dln(p)/dlnr, where (p) is obtained by 
stacking halos in the highest and lowest Zip quartiles. 
We also show the slope expected from NFW as dotted 
curves. The slopes deviate substantially from NFW for 
most radii. 

The effect of formation time on the halo concentration 
can also be seen in the left panel of Figure [8] We mark 



the location of r s as vertical dashed lines. For early- 
forming halos, r s is exactly the radius where the slope 
equals -2. In contrast, for late-forming halos, r s is smaller 
than the radius where the slope equals -2. This indicates 
that the NFW model does not provide an adequate fit 
to their profiles, leading to the correlation between Zi/ 2 
and Anfw seen in Figure [4] In addition, we note that a 
lower r around 0.2i? V j r corresponds to a larger r s , and 
thus a smaller NFW concentration. 
In the notation of the Einasto profile, 



dlnT 
dlnr 



r-2 



^Einasto ■ 



(18) 
(19) 



For early-forming halos, T is close to a power law; there- 
fore, the Einasto profile provides a fit better than NFW^] 
However, for late-forming halos, neither model provides 
a good fit, and the logarithmic slope of T has a sudden 
increase around 0.3 i? v j r . This "kink" in T can be ex- 
plained by the presence of subhalos. In the right panel 
of Figure |HJ we present the slope with the massive sub- 
halos (v max > 200 km/s) removed. After removing these 
subhalos, the kink in the red curve disappears, and the 
difference between the blue and red curves is significantly 
reduced^] We also note that in the presence of a massive 
subhalo, the slope can be shallower near the locus of 
the subhalo and steeper at larger radii. Therefore, the 
presence of subhalos can explain the kink in the slope in 
Figure [HJ 

4.2.3. Evolution of concentration 

Figure [9] shows the mean concentration evolution of 
our main nalos. The red/blue curve corresponds to the 
low/high zi/2 quartile. To emphasize the late-time phe- 
nomenon, we plot logarithmic scale factor as the x-axis. 
For a < 0.6, the averages of both populations remain 
relatively constant and have similar values, reflecting the 
split in Zi^ 2 . For a > 0.6, the concentration of early- 
forming halos increases steadily with time and the scatter 
decreases; this steady increase in concentration is presum- 
ably related to the lack of significant mass growth and 
the resulting higher degree of relaxedness. On the other 
hand, the mean concentration of late-forming halos does 
not increase much and still shows large scatter. 

The right panel of Figure [9] shows the evolution of r s in 
physical /i -1 Mpc. When the mean concentration of early- 
forming halos increases, the mean r s value remains nearly 
constant, suggesting that the increase in concentration 
could be mainly driven by the pseudo-evolution of i? V j r 
related to the time dependence of A virj o cr j t . To confirm 
this, we add an inset showing the mean mass evolution 

9 We calculate the residual for Figure [8j ^nodel = 
J2iLi |rtrue(r») - r modi , 1 (T-i)|/(r t ruc>. For early-forming halos, 

0.65. For late-forming halos. 



2.9. For all halos, 



1.59 



A NFW = 12 and A Eina.to 

A nfw = 3-0 and Ag lnasto = 
and A£. . = 1.02. 

Einasto 

Our subhalo removal procedure is based on the subhalo-particlc 
assignment based on ROCKSTAR; different halo finders tend to have 
different subhalo-particle assignment criteria (e.g., see |Onions et ah] 
2012). Since our purpose is to understand the trend with and 
without the presence of the subhalos, we do not expect that detailed 
particle assignment will impact the trend discussed here. 
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Figure 8. Left: The absolute value of the logarithmic slope of the density profile, F The red/blue curve corresponds to the mean value of 
the low/high quartile. We also mark the r s and slope expected from NFW. Except for the early-forming halos at around r a , F deviates 
substantially from NFW. For late-forming halos, T < 2 at r s (shallower than NFW expectation). For early-forming halos, V approximately 
follows a power law and can be well described by an Einasto profile. For late-forming halos, T appears to be a broken power law and 
suddenly increases around 0.3 -R v ir- Right: Halos with their massive subhalos (v max > 200 km/s) removed. For late-forming halos, the kink 
around 0.3-R v j r disappears with the removal of subhalos, and their differences from the early-forming halos are reduced. This result indicates 
that the differences in T in the left panel are largely driven by massive subhalos. 




Figure 9. Evolution of NFW concentration (left) and scale radius (right). The red/blue curve corresponds to the mean of the halos in 
low/high Z1/2 quartile. For late-forming halos, the mean concentration remains relatively constant for all redshifts. For early-forming halos, 
the mean concentration remains constant at high redshift; however, for a > 0.6, the concentration steadily increases, while r a remains nearly 
constant. In the inset, we compare their late-time mass accretion with the pseudo-evolution of a static halo, indicating that the increase in 
concentration of early-forming halos is largely driven by pseudo-evolution. 
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Figure 10. The mean radial velocity profiles of Rhapsody halos. 
The red/blue curve corresponds to halos in the low/high Z\/2 
quartile, and the band corresponds to the standard deviation of 
the sample. Late-forming halos show a significant outflow between 
0.3 and 1 R v i r on average and have larger scatter. In addition, 
late-forming halos tend to have stronger infall than early-forming 
halos beyond i? v ir- We will demonstrate that the outflow of the 
late-forming halos is related to the coherent motions of dark matter 
particles associated with subhalos. 

of both populations as well as the pseudo-evolution of a 
static halo. As can be seen, for early-forming halos, the 
mass evolution for a > 0.6 is close to the pseudo-evolution, 
indicating that these halos are on average close to static 
and their increase in concentration is indeed driven by 
the pseudo-evolution. We note that the constant r s has 



been shown in lower mass systems (e.g., Bullock et al 



2001 Wechsler et al. 2002), here we show that it also 



happens in relaxed massive halos 

5. THE IMPACT OF FORMATION TIME ON THE 
PHASE-SPACE STRUCTURE 

In this section, we present the impact of formation time 
on the phase-space structure of halos. The phase-space 
structure has been explored for the purpose of seeking 
alternative definitions of halo bound ary and mass (e.g., 
Busha et al.|2005 Cuesta et al.| 2008); understanding the 
connection betw een density and velocit y distr ibution (e.g., 
Hansen fc Moore]|2006j [Navarro et al. 20101; calculating 
possible caustics to aid the s earch for tensing or dark 
matter annihilation signal (e.g., Diemand & Kuhlen 2008). 
In this work, we focus on the radial velocity profile because 
the infall and outflow patterns carry the information of 
formation and merger history and can help us understand 
the build-up of the density profile as well as the state of 
equilibrium. 

Figure [lO] presents the differential radial velocity profile 



i? v 



(20) 



the halos in the highest and lowest quartilespl The late- 
forming halos (red curve) show an average outflow within 
i? v ; r . As will be discussed in the following paragraph, 
these outflows are related to the coherent motions of 
the particles associated with recently merged subhalos. 
When these subhalos pass though the main halo, their 
particles tend to have high outflowing velocities. On the 
other hand, early-forming halos tend to show more regular 
infall patterns. 

The presence of a significant fraction of particles ex- 
ceeding the virial velocity of the main halo is consistent 
with the detection of the so-called "back-splash" popu la- 
tion of halos (e.g., |Gill et al.|2005| [Ludlow et al.|[2009| in 
the outskirts of massive halos. These halos have left the 
virial radius after their initial accretion onto the halo and 
have experienced tidal stripping much like their siblings 
remained inside the virial radius. 

Figure [TT] compares the phase-pace structures for early- 
forming and late-forming halos. We stack the halos in 
the high (left panel) and low (right panel) Z\/% quartiles 
and present V r vs. r (both normalized using the virial 
units). These two populations show different features in 
the phase-space diagram. The late-forming halos again 
show an average outflow between 0.3 and 1 i? V ir- In 
addition, in this radial range, the shapes of the contours 
are different. For early-forming halos, the contours are 
close to parallel to the cc-axis; for late-forming halos, the 
contours are curved upward. 

We find that the outflow in late-forming halos, like the 
kink in T, is related to the pres ence of massive subhalos. 
In the bottom panels of Figure [TTJ we remove subhalos 
with v max > 200 km/s, and the outflow is reduced. 

Although formation history does not st rong ly impact 
the velocity ellipsoid (as will be shown in [6.21, it clearly 
impacts the radial velocity profile. We also note that 
the outflows shown here are more likely to be found in 
dynamically young systems, including cluster-size systems 
with late formation time (like our sample) or galactic 
systems at high redshift. 

6. SHAPES AND ALIGNMENTS 

In this section, we present the shapes for the spatial 
distributions and velocities of dark matter particles of the 
main halos in Rhapsody, as well as the alignments of 
their orientations. Our motivation is again related to the 
formation history: The accretion of matter onto massive 
halos is presumably correlated with the surrounding large- 
scale structure, and anisotropic accretion onto the halo 
can leave an imprint in both the shape of the halo and 
the motion of matter in its interior. 

6.1. Halo shapes 

N-body simulations have shown that dark matter halos 
are generally not spherical ob jects but have significant 
ellipticity and tri axiality (e.g. 
k, Evrardl [20051 jAllgood et a 



Jing fc Suto||2002| |Kasun| 
\ \M&l IBett et 



Hayashi et al.||2007). Calibrations of halo shapes from 



simulations directly im pact the accuracy of weak lensing 
mass calibration (e.g. , Corless & King 2007 Bett 2012 
Schneider et al.|[2012 ) and can be used to constrain the 



of the dark matter particles that are bound to the main 
halos. We again split halos by their Z\/2 and present 



11 We note that for each halo, (V r ) is itself the mean of the 
particles in a radial bin, and the stacking process is the average 
over the (V r ) of individual halos. 
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Low Z]y 2 quartile 



0.1 



0.5 
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High Zj/2 quartile; subhalo-removed 




Low zj/2 quartile; subhalo-removed 




Figure 11. The impact of formation time on the phase-space structure of halos. Upper/Lower: Including/excluding subhalos. Left/Right: 
The high/low Z1/2 quartile. Comparing the two upper panels, we can see that the late-forming halos show stronger outflow between 0.3 and 
1 -Rvin and the curvy shapes of the contours in this region correspond to the kink in the slope of density profile at the same scale shown in 
Figure [8] Comparing the two right panels, we can see that the outflow in the subhalo-removed systems is reduced. 



self-interaction cross-s ection of dark matter particles (e.g., 
Miralda-Escudi|p002| . 

The shape parameters are typically defined through the 
mass distribution tensor with respect to the halo center 



h 



(nrj) 



(21) 



where is the i component of the position vector r of a 
dark matter particle with respect to the halo center, and 
the average (•) is over all dark matter particles within 
i? v ir- Since all particles within R v - lr are of the same mass, 
no weighting by mass is needed. The eigenvalues of 
are sorted as Ai > A 2 > A3, and the shape parameters 
are defined as a = y/Xi, b = y/M, c = y/Xs. We present 



the dimensionless ratios b/a and c/a. In addition, the 
triaxiality parameter is defined as 



T = 



(22) 



We list these halo shape properties in Table [2] and 
present the distribution of c/a in Figure [4] The shape 
parameter c/a is only weakly correlated with formation 
history. The correlation is even weaker for b/a and T. 

It has been shown that the halo s hape depends on the 
radiu s at which it is measured (e.g., Bailin & Steinmetz 
2005 1 . We also measure the halo shape at -K500C and find 



that b/a, c/a, and T measured at R§oo c are correlated 
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with those measured at i? v ir, with correlation coefficients 
0.61, 0.63, and 0.60, respectively. The major axes mea- 
sured with i?5ooc and i? v ir have a median an gle of 21°, 
which agrees with the value (~ 20°) reported by Schneider 
| et al.| (2012) for the Millennium Simulations 



sured along different lines of sight can be calculated as 



|Allgood et al 



(2006) used a smaller radius 0.3i? v j r and 
lod to calculate the reduced inertia tensor 
2 to suppress the influence of the larger 



an iterative met 
(weighted by r~ 
radii). They found that c/a is correlated with formation 
time, and the correlation is weaker for higher mass. In 
our measurement with unreduced inertia tensor at i? v j r 
and -R5000 the correlation between c/a a nd Zi/ 2 is not 
strong (0.27 and 0.29 respectively). Since Allgood et al. 
(2006) did not state the correlation coefficient and did 
not have statistics in our mass regime, we cannot make 
a direct comparison but note that the different mass 
regime and measuremen t methods could impact these 
correlations. Bett (2012) recently showed that the shape 
measured from the iterative reduced tensor is similar to 
and only slightly more spherical than those obtained with 
t he simple me t hod a pplied here. 



Shaw et al. 



(120061) showed that the measured shapes 
depend somewhat on the state of relaxedness. For our 
sample of halos, we also find a weak correlation between 
c/a and the goodness-of-fit proxy A NFW . While the corre- 
lation is weak, it is however interesting to note that none 
of our halos is close to spherical (high c/a) and has a high 
Anfw at the same time. Given such a correlation, one 
might thus wonder whether the deviations from the NFW 
profile — with a strong depe ndence on the formation time 
Zx/2, as discussed in Section |4 . 2 . 3|— are driven by system- 



atic deviations from sphericity of the entire virialized halo 
or by the anisotropic distribution of the most massive 
subhalos. To decide this, we repeat the measurements of 
the shape parameters after removing all subhalos with 
Vmax > 200 km/s. We find that the correlation coeffi- 
cient between Anfw an d c/a at i? v j r drops from —0.20 
to a mere 0.01 when subhalos are removed. Similarly, 
at i?5ooc we observe a drop of the correlation coefficient 
from —0.29 to 0.05. This is strong evidence that the 
deviations from NFW at variance with formation time 
are predominantly driven by recently accreted subhalos. 
The correlation between shape and formation time then 
is a simple consequence of the anisotropic distribution of 
these subhalos. 



6.2. 



White et al. 



Velocity ellipsoid 
have 



(2010) have demonstrated that the 
anisotropic motion ot subhalos in clusters introduces sig- 
nificant scatter in the velocity dispersions measured along 
different lines of sight. Here we follow the same procedure 
to measure the properties of the velocity ellipsoid of the 
dark matter particles in the Rhapsody sample. 

Analogous to the mass distribution tensor, the velocity 
ellipsoid is defined as 

4 - M ■ (23) 

Sorting the eigenvalues of the velocity ellipsoid as Ai > 
A2 > A3, one can again define shape parameters a^"' = 
y/Xi, M") = V^2i c — "v/Als, and dimensionless ratios 
frW/aW and c^/a^ to describe the velocity ellipsoid. 
The mean and scatter of the velocity dispersions mea- 



(^L) = J(Ai+A 2 + A 3 ) , 

2 \2. ' 



(24) 



{So-LY = ^(K + K + H- AiA 2 - A 2 A 3 - A3AO . (25) 

We list these parameters in Table [2] In Figure |4j we 
show the distribution of c^/a^ v \ which, like c/a, is only 
weakly correlated with formation time. The correlation 
is similar for and Serf . 

We also measure the velocity ellipsoid based on R^onc- 
In this case, the correlation is stronger than halo shape: 



h (v)/ 



and Saf Q 



measured at -R500C are corre- 



lated with those values measured at i? v ir with correlation 
coefficients 0.80, 0.86, and 0.86, respectively. The ma- 
jor axes of the velocity ellipsoid measured with i? v j r and 
i?500c have a median angle of 14° . 

As expected, the ellipsoid of dark matter distribution 
and velocity ellipsoid are correlated. While b/a and 
Jo/v) have a correlation of 0.35, we find a stronger 
correlation between c/a and /a^ of magnitude 0.53. 
We note that the velocity ellipsoid is in general more 
spherical than velocity ellipsoid, this difference is related 
to the fact that the shape of the gravitational potential is 
in general more spheri cal than density distribu tion. This 
trend was discussed in |Kasun fc Evrard| ( 2005 ) , who used 
the Hubble volume simulation and denned halos based 
on A 2 ooc- Although our simulation uses a different mass 
definition and is in a different mass regime from theirs, 
our axis ratios for spatial distribution and velocity are in 
agreement within 0.01 with theirs. In addition, the major 
axes of the spatial distribution and velocity ellipsoids 
have a median a ngle of 27°, which is slig htly larger than 
the value 22° inlKasun & Evrardl (12005} . 



We do not find significant correlation between any of 
these shape parameters with the environmental parame- 
ters we measured (including several halo number overden- 
sity and nearest neighbors for several mass thresholds). 
In Paper II, we will repeat these measurements for the 
spatial distributions and velocities of subhalos and com- 
pare them with the results measured from all dark matter 
particles shown here. 

7. SUMMARY AND DISCUSSION 

We have presented the first results of the Rhapsody 
project, which includes 96 cluster-size halos with mass 
M vir = lO 14 ' 8±o - O5 /i- 1 M , re-simulated from 1 /i~ 3 Gpc 3 
with a resolution equivalent to 8192 3 particles in this 
volume. In addition to achieving high resolution and 
large statistics simultaneously, Rhapsody is unique in 
its well-resolved subhalos and wide span of evolution 
history. Rhapsody also implemented the state-of-the-art 
algorithms for initial conditions (Music), halo finding 
(Rockstar), and merger trees. 

Our findings are summarized as follows: 

1. Properties of the cluster halos: We have sum- 
marized the key properties (including various mass 
definitions, formation history proxies, halo concen- 
tration, and shape parameters) of the 96 main halos 
of the Rhapsody sample in Table [2] and shown the 
distributions and mutual covariance of a subset of 
these properties in Figure [ij 
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2. Mass accretion history and merger rate: We 

have investigat ed t he mass accretion history of the 
main halos in S3.1 tracking the most massive pro- 
genitors over 5 decades of mass growth. We have 
shown that an exponential form does not provide 
an adequate fit to such a wide span of t ime, and 
an extra power law is needed. In §3.2| we have 
shown that the differential merger rate follows a 
power law scaling with the merger mass ratio and 
is independent of redshift, in agreement with earlier 
work based on less massive halos or smaller samples. 

3. Density profile: In Sj4j we have quantified in de- 
tail how formation time impacts the halo density 
profile. We have shown that the deviations from the 
NFW model systematically depend on the forma- 
tion history of the halo. Specifically, late-forming 
halos t end to show larger deviations from NFW. In 
{ 4.2.2[ we investigate the slope of the density profile. 
Early-forming halos can be well described by the 
Einasto model. For late-forming halos, the slope 
appears to be a broken power law with a sudden 
transition around 0.3 R V i r . This kink in the slope 
profile is related to the abundant massive subhalos 
i n thes e late-forming halos. We have also shown in 
§4. 2. 3| that, for early-forming halos, the evolution 
of concentration between z ~ 1 and is consistent 
with pseudo-evolution. 

4. Halo shape and alignment: As discussed in Sj6j 
the shape and velocity ellipsoid of our sample is only 
moderately correlated with the formation history. 
We present the correlation and alignment of shapes 
and velocity ellipsoids measured at i? v ; r and i?5oo c - 
Our results indicate that anisotropy of the entire 
virialized halo alone cannot explain the deviations 
of the density profiles from the NFW model and its 
dependence on formation time. 

5. Phase-space structure and formation time: 

We have provided evidence in £j5]for a connection 
between the density profile, massive subhalos, and 
the formation history of the halos. When investigat- 
ing the r-V r phase space, we find that late-forming 
halos show evidence for large fluctuations in the 
radial velocities of their dark matter particles, with 
a considerable fraction in excess of the virial veloci- 
ties, leading to an outflow of mass from these cluster 
halos. We associate fluctuations in the radial ve- 
locities, outflows, as well as the kink in the density 
profile slope with recently accreted subhalos. On the 
contrary, early-forming halos show a substantially 
more regular motion of dark matter, as expected 
if they are dynamically more relaxed. Our results 
thus indicate that the recently accreted subhalos 
drive the observed deviation from the NFW profile. 
Due to their anisotropic distributions, these sub- 
halos also lead to a correlation between formation 
time and halo shape. 



In Paper II, we analyze in detail the properties of the 
subhalo population of the cluster halos. In addition to 
the impact of formation time on subhalo properties, we 
investigate the intertwining effects from subhalo selection, 
stripping, and resolution of simulations. 

The Rhapsody simulation suite provides valuable in- 
formation for other aspects of cluster cosmology. For 
example, the cluster-size halos in Rhapsody can be fur- 
ther used to study the covariances between mass tracers, 
for example, galaxy content, dynamics of galaxies, weak 
gravitational lensing, X-ray, and the SZ effect. The forma- 
tion history and environment of clusters can potentially 
impact these mass proxies systematically, either by al- 
tering the intrinsic properties of clusters or by affecting 
the observed properties through the line-of-sight projec- 
tion. As current multi-wavelength surveys combine these 
different observables for cluster mass calibration, it is 
imperative to understand the covariances between these 
observables. 

Our re-simulation technique can also be applied to study 
the "pink elephant" clusters, which refer to a handful of 



Jee et al. 


2009 


Foley et al. 


2011) 



stimulated a great amount of discussion about whether 
they pose a c hallenge to the current ACDM paradigm 
of cosmology ( |Mortonson et al.||2011| |Hoyle et aL]|2011[ ). 
To interpret the cosmological implications of these clus- 
ters correctly, it is important to understand their mass 
calibration. An extension of the current Rhapsody sam- 
ple that includes a statistical sample for these massive 
clusters at high redshift will improve our understanding 
of these massive clusters. Understanding the covariances 
between different observable quantities and the potential 
biases in the mass measurements of these clusters can 
help us disentangle the astrophysical and cosmological 
implications of these clusters. 



This work was supported by the U.S. Department of 
Energy under contract numbers DE-AC02-76SF00515 
and DE-FG02-95ER40899, and by Stanford University 
through a Gabilan Stanford Graduate Fellowship to HW 
and a Terman Fellowship to RHW. Additional support 
was provided by SLAC-LDRD-0030-12. This work uses 
data from one of the LasDamas simulations; mock 
galaxy catalogs from these simulations are available at 
http:/ /lss. phy.vanderbilt.edu/lasdamas/ . We are grateful 
to our LasDamas collaborators, and especially Michael 
Busha, who ran the Carmen box that was used for re- 
simulation; we also thank Michael for extensive helpful 
discussions. We thank Ralf Kaehler for assistance with 
visualizations of the simulations. We gratefully acknowl- 
edge the support of Stuart Marshall, Amedeo Perazzo, 
and the SLAC computational team, as well as the compu- 
tational resources at SLAC. We also thank Gus Evrard, 
Eduardo Rozo, Matt Becker, Andrey Kravtsov, Sarah 
Hansen, Anja von der Linden, Douglas Applegate, and 
Will Dawson for helpful discussions and comments. 



18 



WU ET AL. 



6 8 10 12 



1 °- 6 

s 




(II). 211. 10. 81. 01. 2 



exp 

exp+PL (2-para) 
exp+PL (1-para) 
exp+PL (1-para, optimal) 



A (RMS of Residual) 



1 




.0 0.8 1.0 1.5 2.(1 2.5 3.0 



Figure 12. Left: The mass accretion rate of the RHAPSODY halos, — dlnM v i r /dz, as a function of z, averaged over every 3 output 
redshifts to reduce the noise. Middle: Comparison of different fitting forms of mass accretion history: exponential model (Eq.j3l with one 
free parameter a; the exponential-plus-power law model with two free parameters /3 and 7 (Eq.|6j or one free parameter (Eq. [9J. The 
one-parameter model using an optimal relation between /3 and 7 (the purple curve; also see the inset) provides a compelling description of 
the mass accretion history. Right: The cumulative distribution of several proxies for formation time: 21/2 (the exact half-mass redshift); 
z a = ln2/a; z y _p = ln2/(7 — zp^ (half-mass redshift obtained by solving Eq.[6|. These different formation time proxies probe somewhat 
different epochs in a halo's history. 



APPENDIX 
FITTING THE FORMATION HISTORY 
For each mass accretion history model described in §3.1| we minimize the target function 



A 2 

model 
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-Y 

i=l 
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logio 



M, 
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We note that this function differs from what was used in McBride et al. ( 2009 1 



A 



2 
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N 2s 

i=l 



M( Zi )/M - M modc i(z t )/M ] 
M(zi)/M 



(Al) 



(A2) 



Because our mass accretion history spans approximately 5 orders of magnitude in mass and starts from redshift 12, 
weighting by M(zi) /Mq si gnifi cantly underweights high redshift outputs. 
The left panel of Figure 12 shows the mass growth rate defined as 



dlnM v 
dz 



(A3) 



as a function of z. In [3.1 we have discussed the deviation of mass accretion history from a power law, and this figure 
clearly shows that the mass accretion rate is not constant and presents a large amount of scat ter. 

The middle panel presents the cumulative distribution of RMS residuals, A mo d c i (Eq. |Al[ ), for the fits to various 
models of accretion history. The red curve corresponds to the exponential model, which lias the largest residuals; 
the green curve corresponds to the exponential-plus-power law model with two free parameters (/3, 7), which has the 
smallest residuals. In between these two models, we compare two 1-parameter models: the blue curve corresponds to 
the model that uses the linear fit between (3 and 7 (/3 = 4.I67 — 4.00) to eliminate one parameter; the purple curve 
corresponds to an optimization of the relation between /3 and 7 to minimize the overall residuals, obtained with several 
iterations. The optimal relation is given by 

(3 = 5.277 - 4.61 (A4) 

This 1-parameter model performs almost equally well as the 2-parameter model does. 

The right panel of Figure |T2| shows the cumulative distribution functions for several halo formation time proxies: 

• Z1/2, the redshift when the halo first reaches half of its final mass. 

• z a = ln2/a, the redshift at which M(z a ) — Mo/2 for the exponential fit. 

• z 7 _^ = In 2/(7 — /3), analogous to z a , where /3 and 7 come from a fit to the exponential-plus-power law model 
(Eq. [6|). We note that z 7 -/3 is equivalent to the value of z a measured with only the low redshift outputs. 

• zpj, obtained by numerically solving M(z^ 7 ) = Mo/2 using the exponential-plus-power law model (Eq. [6]). 
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Figure 13. Left: The Kolmogorov-Smirnov statistics for three models for the density profile: NFW, NFW-like (with a free outer 
slope), and Einasto. The NFW-like and the Einasto models work equally well. Middle and Right: Concentration distributions based on: 
C NFW> c NFW-likei an d cginasto = Rvir/R—2- Middle/Right shows the cumulative distribution of c/lnc and the corresponding best-fit of 
normal/log-normal distribution (dashed curve). Both normal and log-normal can well describe cnjp\v an d CE; nasto , while cnpw— like prefers 
a log-normal distribution. 

None of the formation time proxies obtained from the fitting functions captures the exact half-mass redshift. The 



rank correlation between z a and 



Zl/2 



is 0.55, and that between 

114.8; 



and 



■1/2 



is 0.69. Since we fit for 5 orders of 



magnitude in mass growth (10 to 10 h~ Mq), these functions are not flexible enough to describe the late stage of 
halo evolution. Nevertheless, these different formation time definitions are still useful because they probe different 
epochs in a halo's history; z a tends to be a slightly earlier epoch, while z 1 -$ and z^ 7 tend to be later compared to zi/ 2 - 
Although zp 7 is completely correlated with it corresponds to a slightly earlier redshift. 

FITTING THE DENSITY PROFILE: GOODNESS OF FIT 

For each model described in Q we apply the maximum-likelihood method for fitting the density profiles of individual 
halos. The procedure of finding the maximum-likelihood estimator is to maximize the log-likelihood function over a set 
of parameters, p. The log-likelihood function is defined as 



(Bl) 



where the summation runs over all the N particles, and 



1 

M 



Airr p p (r) 



(B2) 



so that J v p (r)dr = 1. This approach is consistent with the radially- binned fitting method with a large number of 
particles and is more stab le th an using a fixed number of bins when the halo has fewer particles. 
The left panel of Figure 13 presents the Kolmogorov-Smirnov statistics for these three models, where 



A 



model 



max|M(< r) -M mode i(< r)\ 

Mvir 



(B3) 



The two-parameter models, the NFW-like and the Einasto models, work almost equally well and are significantly better 
than the single-parameter NFW model. 

In the middle/right panel, we show the cumulative distribution for c/lnc for the three different models stated 
above. We also show the corresponding best fit normal/log- normal distributions and list the p-value based on a 
Kolmogorov-Smirnov test for goodness of fit. For the NFW and the Einasto models, both normal and log-normal 
distribution provide acceptable descriptions. For the NFW-like profile, a log-normal distribution provides a slightly 
better description. 
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